f=@(x)x-1;
t = fzero(f,0);


plot_RAS()
function plot_RAS()
    theta = linspace(0, 2*pi, 1000);
    z = exp(1i * theta);
    figure;

    % === Adams-Bashforth (order 1–3) ===
    subplot(2,2,1);
    hold on;
    for p = 1:3
        [rho, sigma] = ab_coeffs(p);
        kappa = rho(z) ./ sigma(z);
        plot(real(kappa), imag(kappa), 'DisplayName', ['AB p=' num2str(p)]);
    end
    title('Adams-Bashforth RAS');
    xlabel('Re(\kappa)'), ylabel('Im(\kappa)');
    axis equal, grid on, legend;

    % === Adams-Bashforth (order 4) ===
    subplot(2,2,2);
    hold on;
    for p = 4
        [rho, sigma] = ab_coeffs(p);
        kappa = rho(z) ./ sigma(z);
         % Create boundary of stability region
        pi_poly = @(zeta) rho(zeta) - kappa * sigma(zeta);
        plot(real(kappa), imag(kappa), 'DisplayName', ['AB p=' num2str(p)]);
    end
    title('Adams-Bashforth RAS');
    xlabel('Re(\kappa)'), ylabel('Im(\kappa)');
    axis equal, grid on, legend;
    

    % === Adams-Moulton (order 3–5) ===
    subplot(2,2,3);
    hold on;
    for p = 3:5
        [rho, sigma] = am_coeffs(p);
        kappa = rho(z) ./ sigma(z);
        plot(real(kappa), imag(kappa), 'DisplayName', ['AM p=' num2str(p)]);
    end
    title('Adams-Moulton RAS');
    xlabel('Re(\kappa)'), ylabel('Im(\kappa)');
    axis equal, grid on, legend;

    % === Backward Differentiation Formulas (BDF 1–4) ===
    subplot(2,2,4);
    hold on;
    for p = 1:4
        [rho, sigma] = bdf_coeffs(p);
        kappa = rho(z) ./ sigma(z);
        plot(real(kappa), imag(kappa), 'DisplayName', ['BDF p=' num2str(p)]);
    end
    title('Backward Differentiation RAS');
    xlabel('Re(\kappa)'), ylabel('Im(\kappa)');
    axis equal, grid on, legend;
end

function [rho, sigma] = ab_coeffs(p)
    coeffs = {
        [1 -1],          [0 1];          % p = 1
        [1 -1 0],        [0 3/2 -1/2];   % p = 2
        [1 -1 0 0],      [0 23/12 -16/12 5/12]; % p = 3
        [1 -1 0 0 0],    [0 55/24 -59/24 37/24 -9/24]; % p = 4
    };
    a = coeffs{p,1}; b = coeffs{p,2};
    rho = @(z) polyval(flip(a), z);
    sigma = @(z) polyval(flip(b), z);
end

function [rho, sigma] = am_coeffs(p)
    coeffs = {
        [1 -1],              [1/2 1/2];            % p = 2
        [1 -1 0],            [5/12 8/12 -1/12];    % p = 3
        [1 -1 0 0],          [9/24 19/24 -5/24 1/24]; % p = 4
        [1 -1 0 0 0],        [251/720 646/720 -264/720 106/720 -19/720]; % p = 5
    };
    a = coeffs{p-1,1}; b = coeffs{p-1,2};
    rho = @(z) polyval(flip(a), z);
    sigma = @(z) polyval(flip(b), z);
end

function [rho, sigma] = bdf_coeffs(p)
    coeffs = {
        [1 -1],              [1];            % p = 1
        [1 -4/3 1/3],        [2/3];          % p = 2
        [1 -18/11 9/11 -2/11], [6/11];       % p = 3
        [1 -48/25 36/25 -16/25 3/25], [12/25]; % p = 4
    };
    a = coeffs{p,1}; b = coeffs{p,2};
    rho = @(z) polyval(flip(a), z);
    sigma = @(z) polyval(flip(b), z);
end




